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We build upon previous work that used coherent states as a measurement of the local phase 
space and extended the flux operator by adapting the Husimi projection to produce a vector field 
called the Husimi map. In this article, we extend its definition from continuous systems to lattices. 
This requires making several adjustments to incorporate effects such as group velocity and multiple 
bands. Several phenomena which uniquely occur in lattice systems, like group-velocity warping and 
internal Bragg diffraction, are explained and demonstrated using Husimi maps. We also show that 
scattering points between bands and valleys can be identified in the divergence of the Husimi map. 



INTRODUCTION 



In Mason et aZ. [lj , we introduced a new interpretation 
of the probability flux operator 
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by expressing its eigenstates as the limit of measure- 
ments by infinitesimally small coherent states. Our ap- 
proach yields a new perspective on flux measurements 
and provides a novel tool for visualizing wavefunctions 
which parallel the probability flux map. Because they 
are based on the Husimi projection technique [2 J, these 
visualizations are called "Husimi maps". Husimi maps 
improve the understanding of the semiclassical paths un- 
derlying the quantum wavefunctions and can be of use 
even for systems where the traditional flux has little suc- 
cess (i.e., when it is either zero or strongly misleading). 
Later work[3j further developed the numerical framework 
of the Husimi map and applied it to a wider variety of 
systems by incorporating local potentials and examining 
flux through open systems. 

This article expands the Husimi map technique 
from continuous, free-particle systems like the two- 
dimensional electron gas (2DEG) to crystalline systems 
like graphene. While the extended wavefunction of an 
electron in a crystal is continuous, the potential imposed 
by the nuclei can be modeled by replacing the contin- 
uum with localized wavefunctions centered at individual 
tight-binding lattice sites. These individual wavefunc- 
tions combine to form a model of the entire wavefunc- 
tion, which now defines their envelope function. In this 
model, Eq. [T] describes not the probability flow at an in- 
finitesimal point, but the flow of probability in and out 
of the localized wavefunction at a single site. 

Lattice systems can behave very differently from con- 
tinuous systems. For instance, the orientation of the 
group velocity vector, which dictates classical dynamics, 
can strongly diverge from the wave vector, which was the 
initial foundation of the Husimi projection. In fact, the 
group-velocity space can be so strongly restricted that 
classical trajectories are only permitted along certain di- 
rections, dramatically affecting the dynamics of states 
that inhabit lattice systems. When these trajectories hit 



a boundary, internal Bragg diffraction can produce addi- 
tional nonclassical ray reflections. 

Here we explore two-dimensional square and honey- 
comb lattices; extension to three-dimensional systems is 
straightforward. Honeycombs induce an additional phe- 
nomenon: the presence of multiple bands and valleys, by 
which different quasiparticles can propagate and inter- 
fere. While the flux operator is unable to reflect any of 
these behaviors, in this article we show that with proper 
modifications, the Husimi projection can handle them 
with ease. 

This paper is organized as follows: In Section |II A| 
we provide the definition of the Husimi projection for 
continuous system and then modify the Husimi projec- 
tion in Section [II C| to represent the group velocity and 
multiple bands. In Section |III A[ we apply the Husimi 
projection to square lattices near the band center where 
group- velocity effects are strongest, and in Section [TlIB[ 
we examine the graphene honeycomb lattice. Finally, 
in Section [HI D[ we provide an interpretation of unusual 
boundary reflections found in Sections III A and IIIB| by 
demonstrating and measuring internal Bragg diffraction. 



II. METHOD 

A. Definition of the Husimi Projection 

Building off work in Husimi [2 J and Mason et al [3], we 
define the Husimi function as a measurement between 
a wavefunction ip ({r^}) and a coherent state |ro,ko,<r), 
which minimizes joint uncertainty in spatial and momen- 
tum coordinates. For lattice systems, the wavefunction 
represents the probability amplitude multiplier of local- 
ized wavefunctions indexed at discrete lattice sites, which 
are associated with discrete positions in the set {r^}. The 
coherent state is also an envelope function over localized 
wavefunctions, defined by the Gaussian function 

e -(ri-ro) 2 /4cr 2 +ikoTi 

centered around i*o and ko. The parameter a defines the 
spatial spread of the coherent state and the uncertainties 
in space and momentum according to the well-known re- 
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As a result, there is a trade-off for any value of a selected: 
for small a, spatial resolution is improved at the expense 
of resolution in momentum space, and vice versa for large 
a. 

We can explicitly write out the projection of the wave- 
function onto the coherent state as 
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where d is the number of dimensions in the system. The 
Husimi function is then defined as 



Hu(r ,k ,cr;^({ r ;})) = |(^| r , k , a) 
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If we weight the Husimi function by the central wavevec- 
tor ko, we obtain the Husimi vector. When momentum 
space is explored at a point by many Husimi vectors, the 
result is the full Husimi projection. 

If all the Husimi vectors at a point are summed, the 
Husimi function can be used to construct and generalize 
the flux operator, resulting in the vector- valued function 
Hu (r , a; ^ ({r*})) equal to 

Hu(r ,cr;^({ri})) = j r o, k , cr)\ 2 k dk . (5) 

Earlier work has shown that for ak <C 1, Eq. [5] is pro- 
portional to the traditional flux vector expectation value 
[3J. For lattice systems, the traditional flux becomes a 
finite- difference approximation defined by the lattice. 



B. The Hamiltonian 

This paper examines Hamiltonians using the nearest- 
neighbor tight-binding approximation 
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where aj" is the creation operator at orbital site i and we 
sum over the set (ij) of nearest neighbors. The quantity 
€i is the on-site energy and t is the hopping energy scale. 
For the square lattices, we set e = —At. For systems 
at energies E < 0.5t, the tight-binding Hamiltonian is a 
close approximation to the effective mass envelope func- 
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tion Hamiltonian H = — ^+U{y) where t = h 2 /(2m*a 2 ) 
and a is the mesh lattice spacing. For the honeycomb 
lattice, parameters are set to the common values in the 
literature for graphene: e = and t = 2.7eV[U[5]. Eigen- 
states of closed stadium billiard systems are computed 
using the standard sparse matrix eigensolvers. 
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Figure 1: The two-dimensional dispersion relation for the 
square lattice demonstrates strong group- velocity warping at 
the band center (E — At). The dispersion relations for E — 
0.9t, and 7. It (dashed white lines) are nearly circular, while 
the relation near the band edge at E — 3.9t (solid white lines) 
shows strong warping. 



C. Group Velocity 

In the original introduction of the Husimi map[3j, each 
Husimi function is weighted by the wavevector of the co- 
herent state to produce a visual guide to the classical dy- 
namics of the system. Summing all the vectors equates 
to the flux operator (Eq. |5| when the coherent states are 
sufficiently small. This equivalence holds in lattices, how- 
ever, the direction and magnitude of the group velocity 
Vk^ (k) can strongly diverge from the wavevector. Since 
a coherent state, which is now defined as an envelope 
function over localized wavefunctions, follows the group- 
velocity vector instead of its wavevector, it is necessary to 
weight the Husimi function by group- velocity vectors to 
indicate the classical dynamics. As a result, the Husimi 
projection indicates the classical flow of quasiparticles, 
in contrast to the flux operator (Eq. ??), which instead 
indicates the flow of probability. 

At low energies, the square lattice closely approxi- 
mates a free-particle continuous system so that this mod- 
ification is minimal. At higher energies, however, the 
mapping from the wavevector to group velocity can be 
strongly constricted. For example, at energies near the 
band center of E = 4t, there are only four directions 
available to the group velocity in the square lattice, as 
shown in the solid white contour in Fig. [T]at E = 3.9t. 

To visualize this effect, we show group- velocity Husimi 
projections at three representative energies in Fig. [2] 
for the square lattice. Thirty- two equally-spaced angles 
along a circle are chosen to represent the local momen- 
tum space. Wavevectors are chosen with these angles to 
satisfy the dispersion relation for a given energy. 

At energies away from the band center for the square 
lattice, semi-classical trajectories can assume any direc- 
tion, but near the band center they must follow pre- 
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Figure 2: The group- velocity Husimi projection for the square 
lattice is strongly affected by warping at energies near the 
band center (Fig. [I]). Husimi projections are shown for 
the square lattice for the group-velocity representation at 
E = 0.9£(a), 3.0£(b), and 3.9£(c) with relative uncertainties of 
Ak/k = 2 (top) and 50% (middle and bottom). A schematic 
of the dispersion relation contour at each energy is shown at 
the far bottom. The generating wavefunction for each row 
is shown in (d). In the top and middle row the test wavefunc- 
tion is a cosine wave pointing along the 45° diagonal, and in 
the bottom along the 0° horizontal. 



ferred directions determined by the group- velocity warp- 
ing. However, the manner in which they do so may differ. 
This can be seen in Fig. [2] which examines two cosine- 
wave states with different wavevectors. As the energy 
of the system increases from left-to-right, group- velocity 
warping draws Husimi vectors, and the classical paths, 
towards four preferred directions. When the generating 
wavevector points along one of these directions, group- 
velocity warping merely sharpens the profile. When the 
generating wavevector points in between the preferred di- 
rections, as in the bottom row of Fig. |2j the classical tra- 
jectories are more strongly dependent upon the system 
energy. 

Any expectation value over a wavefunction must be 
evaluated, usually by an integral, over a complete ba- 
sis. As a result, any expectation value derived from 
the Husimi projection must be first computed from the 
wavevector basis; modifications to account for group 
velocity are determined afterwards. For instance, in 
the Multi-Modal Algorithm, which approximates the 
full Husimi projection by a subset of local plane waves 
(see Mason et al.[3\ for more details), each approxima- 
tion is achieved by computing the dot product between 
the Husimi projection and template projections in k- 
space. Because of group-velocity warping, the result- 
ing wavevectors no longer indicate classical flow. To ad- 
dress this problem, the resulting dominant wavevectors 




Figure 3: Like the square lattice in Fig.[T] the two-dimensional 
dispersion relation for the honeycomb lattice demonstrates 
strong group- velocity warping at energies away from the Dirac 
point. The dispersion relation for E = 0.5£ (dashed white 
lines) is nearly circular, while the relation at E — 0.98£ (solid 
white lines) shows strong warping. The K and K' valleys are 
indicated. 



are then mapped onto group velocity by taking the local 
derivative of the dispersion relation. 



D. Band Structure 

The number of bands for a lattice system is equal to 
the number of tight-binding orbitals in the unit cell[6j. 
The square lattice has only one unique tight-binding 
orbital and only one band, but due to the warping in 
the band structure, distinct behaviors result at energies 
above E = 4£, corresponding to the hole pocket (see the 
contour lines in Fig. [T] near the corners of the Brillouin 
zone). However, because the quadratic dispersion rela- 
tions at E = and E = 8£ are separated by energy, a 
semi-classical interpretation of a wavefunction is always 
constrained to one relation or the other. 

In the honeycomb lattice, however, there are two 
unique orbitals in the lattice structure, yielding two 
bands that touch at the Dirac point at E = 0£. But 
more interestingly, the band structure warps each band 
to produce the inequivalent K and K' valleys at the 
Dirac point, which are indicated in Fig. [3] Unlike the 
square lattice, these two valleys co-exist in the energy 
range —t<E<t. 

These valleys exhibit a linear dispersion relation near 
the Dirac point |4J. At energies away from the Dirac 
point, the two valleys undergo group- velocity warping 
that emphasizes three directions, which is referred to as 
"trigonal warping" |5J. The effects of trigonal warping can 
be significant even at energies as low as 0.2£. 

The Husimi map can assist visualizations of intervalley 
scattering, the scattering of quasiparticles between the 
two valleys which are part of the same band. To resolve 
the two valleys, it is simply necessary that the uncer- 
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Figure 4: The full wavevector Husimi map (left), multi- 
modal analysis (middle) and wavefunction (right) for two sta- 
dium eigenstates at energies E\ = 1.496£ (a) and E 2 = 3.982£ 
(b) (a schematic of the dispersion relation contour at each 
energy is shown in the insets). The uncertainty for each pro- 
jection is set to Ak/k = 10%, and the spread of the coherent 
state is indicated by double arrows on the right. Angular 
deflection (Eq. [8} is indicated in blue. Each eigenstate has 
similar characteristic wavelengths, but the lower eigenstate is 
sampled with half the linear resolution, causing its energy to 
go up and the group velocity to become more restrictive. 



tainty of the coherent state is small enough in /c-space 
to unambiguously resolve the wavevectors of each valley. 
Because the two valleys of graphene are well-separated 
and only come close at the corners of each triangle in 
Fig. [3j the Husimi projection can clearly resolve the two 
valleys at most energies. More complicated lattices can 
have additional bands, and any automated method for 
calculating Husimi maps for these systems have to take 
their mutual distance in /c-space into account. 



III. RESULTS 



Stadium Billiard Eigenstates of the Square 
Lattice 



In Fig. |4j we examine two closed stadium billiard sys- 
tems with identical geometric parameters. Both systems 
are created using the square-lattice tight-binding model, 
but the lattice constant in Fig. [4Jd is twice as large, 
so the system possesses far fewer sampling points and 
experiences stronger effects from group- velocity warp- 
ing. These systems connect to the lattice-sampled 
Schrodinger equation for a continuous system, which can 
avoid the effects of group-velocity warping by increasing 



the number of sample points in the system (See Fig. [4|. 
However, lattice spacing is not an adjustable parame- 
ter in atomic systems, and group velocity must be given 
careful attention. 

Keeping the characteristic wavelength constant raises 
the energy in any system with a longer lattice con- 
stant. In Fig. [4]3, an eigenstate of the system is 
shown with energy E2 = 3.892£, near the band cen- 
ter. The energy for the system in Fig. [4^i is chosen 
to reflect the same characteristic wavelength, which de- 
pends upon which direction in /c-space is considered. 
Along the fc^-axis, the energy is bounded below by 
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^ = \ and ti = £2 = t an eigenstate is chosen with an 
energy near the average of the bounds at E\ = 1.496t. 

The classical trajectories indicated by the Husimi map 
in the low-energy state in Fig. [4^i point along directions 
oblique to the 45° diagonals and intermingle among other 
paths at other angles. The Husimi map for the higher en- 
ergy eigenstate in Fig. [4Jd instead only indicates classical 
trajectories along the 45° diagonals. Moreover, the tra- 
jectories in the higher-energy system are much clearer 
since they are reinforced by a restricted group-velocity 
space. 

The Husimi map makes it possible to measure "angu- 
lar deflection", which reflects how classical trajectories 
deviate from the straight line in response to the sys- 
tem. Angular deflection thus provides a map of where 
the boundaries or external potentials most strongly ef- 
fect these dynamics, and can be interpreted as a force on 
the particle represented by the wavefunction. 

For lattice systems, the original definition of angular 
deflection provided in Mason et al.[3\ must be modified 
to account for group velocity. It can thus be defined 



Qang. (r;tf) = j Dabs.(r,k; tf) |V k £ (k)| d d fc, 



(7) 



where the quantity Vk-E (k) represents the group- 
velocity vector corresponding to the wavevector k', and 
the integral covers all wavevectors satisfying the disper- 
sion relation. The quantity D a bs.( r 5 k; \&) is defined as 
the Gaussian- weighted absolute divergence of the Husimi 
function for one particular trajectory angle 
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where we sum over the d orthogonal directions each as- 
sociated with unit vector e$. 

Fig. [4] shows angular deflection in blue, concentrated 
on the boundary as expected. Because the resolution of 
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angular deflection is limited by the spread of the coher- 
ent state used for Husimi sampling, its spread into the 
bulk from the boundary exhibits the same Gaussian dis- 
tribution that is used for the test wavepacket. It is worth 
noting that without the proper modifications, angular de- 
flection based on the wavevector shows non-trivial results 
in the bulk of the system even when there are no external 
fields. 

This also suggests that modifications may be in or- 
der for other metrics for lattice systems. For instance, 
by coordinating the boundary divergence with each 
wavevector, one can compute the quantum analog of a 
state's Poincare map[7j. In Birkhoff coordinates [8j |9], 
the angle of impact is mapped against a coordinate 
along the boundary, and both fully quantum[lOl [11] 
and classical [12J variations have become valuable tools 
in quantum chaos. By incorporating group- velocity con- 
siderations, these metrics may be extended to lattice sys- 
tems. 



Stadium Billiard Eigenstates of the Honeycomb 
Lattice 



For the square lattice, time-reversal symmetry is ex- 
pressed in the Husimi projection by the fact that each 
Husimi vector is accompanied by another of equal mag- 
nitude but opposite direction. This causes the flux oper- 
ator and Eq.[5]to return null results. The same is true for 
the honeycomb lattice, except that the range of wavevec- 
tors available at low energies point towards the K and 
K' valleys. 

But when the Husimi vectors are weighted by the 
group- velocity and not the wavevector, a different be- 
havior emerges. In the honeycomb lattice, group- velocity 
doesn't correlate at low energies with k but k — K^. If 
one examines the Husimi projection for each valley in- 
dividually, it is no longer true that each Husimi vector 
is accompanied by its opposite. Rather, each valley is 
the time-reversal symmetric version of the other, allow- 
ing Husimi vectors in each valley to sum to non-trivial 
results. 

Fig. [5^i shows the Husimi map of the K' valley for a 
high-energy eigenstate in part (b) where the strong pull 
towards the three preferred group- velocities is evident. 
In parts (c) and (d) , the multi-modal analysis for the K' 
and K valleys are shown. According to the time-reversal 
symmetric relation, the Husimi map for the K valley is 
the precise inverse of the K' valley. While the classical 
trajectories are evident in the wavefunction (Fig. [5)3), the 
Husimi map identifies their orientation for each valley. 

Because Husimi vectors for each valley no longer sum 
to zero, it is possible to produce divergence in the Husimi 
map for each wavevector. This is identical to angular de- 
flection in Eq. [8] except that the absolute value is not 
taken in the integrand. And like angular deflection, the 
integrand must be weighted by the group- velocity vector, 
or non-trivial results emerge in the bulk of the system. 




Figure 5: The full Husimi map around the K' valley (a), the 
wavefunction(b), the multi-modal analysis for the K' valley (c) 
and for the K valley (d) for a high-energy eigenstate of the 
honeycomb lattice at E = 0.786£. This system is a closed sta- 
dium billiard system with 20270 lattice points. The relative 
uncertainty in all calculations is Ak/k — 20% with the coher- 
ent state spread indicated by the double- arrows. Because of 
time-reversal symmetry, the Husimi maps in (c) and (d) are 
exact inverses of each other. Unlike the square lattice, the 
summing the Husimi vectors for each valley in a honeycomb 
lattice gives non-zero results for a closed system, giving rise to 
non-trivial divergences along the boundary where one valley 
scatters into the other (indicated in green for positive and red 
for negative). 



Summing the divergence for all wavevectors produces 
the total divergence, which appears in green and red in 
Figs. [5^ and|5]i to indicate positive and negative values. 
These points are, in fact, sources and drains for each val- 
ley, and represent the inter- valley scattering points along 
the boundary, whose scattering properties depend on the 
angle of the cut [13]. The results in Fig. [5] suggest that 
each classical trajectory in this wavefunction shares half 
of its existence in one valley, and half in the other. 



C. Group Velocity Warping 

This section expands upon our findings in Figs. [4] and [5] 
by examining Husimi projections in detail. In Fig. [6| we 



(a) , 
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Figure 6: Two stadium eigenstates for the square lattice 
(a) and the honeycomb lattice (b). The wavefunctions 
(left), wavevector Husimi (middle) and group- velocity Husimi 
(right) projections are shown for the points circled in red. Un- 
certainties for both projections are Ak/k = 20%. Not only is 
there more spread to the wavevector projections, these pro- 
jections also indicate markedly different trajectory paths than 
the group- velocity equivalents. Moreover, the group- velocity 
projections are more consistent with the paths indicated by 
the wavefunctions. 

show Husimi projections for the square (a) and honey- 
comb (b) lattices in both wavevector and group- velocity 
representations. As expected, the spread of each Husimi 
projection is dramatically reduced in the group- velocity 
representation, a consequence of group- velocity warping 
and consistent with Fig. [2] Moreover, a close exami- 
nation reveals that Husimi wavevectors can point along 
surprisingly divergent angles from their trajectories, em- 
phasizing the extent to which group- velocity warping es- 
tablishes such states. 

If it is possible to produce similar classical trajectories 
using a wider variety of wavevectors for lattices, then how 
are wavevectors distributed in this wider range? We can 
provide an answer by summing the Husimi projections 
over a range of eigenstates. We find that with a suffi- 
cient range of eigenstates, neither wavevector nor group- 
velocity distributions vary across the bulk of the system, 
except along the boundaries. For the square-lattice bil- 
liards, directions parallel to boundaries are emphasized, 
which is consistent with Dirichlet boundary conditions. 
This occurs even for jagged edges that do not fall along 
a symmetry axis of the underlying lattice. 

Boundaries on the honeycomb lattice emphasize either 
parallel trajectories for intra- valley scatterers (zig-zag) or 
perpendicular trajectories for inter- valley scatters (arm- 
chair) . We find that honeycomb lattice systems are more 
sensitive than the square lattice to the set of states we 
sum over, requiring a larger sum to provide uniform re- 
sults. The distance from the edge where the emphasis 
occurs is a function of the characteristic wavelength; for 
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Figure 7: The distribution of Husimi vectors from the red 
circles in Fig. [6] summed over hundreds of eigenstates near 
E — 3.5t for the square lattice and E — 0.8t for the honey- 
comb lattice, with a coherent wavepacket spread of Ak/k — 
10%. Above, the wavevectors in the square lattice (a), group- 
velocities in the square lattice (b), wavevectors for the K' 
valley in the honeycomb lattice (c), and group- velocities for 
the K' valley in the honeycomb lattice (d). Husimi projec- 
tions tend to emphasize wavevectors away from the preferred 
group- velocity directions (a,c), but not enough to overcome 
that preference in the group- velocity distribution. 



the square lattice, this is k and in the honeycomb lattice 
q = |k — K'|. As a result, for small enough honeycomb 
systems, the emphasis of directions along the boundary 
can persist into the bulk. 

In Fig. [7| we show a representative distribution of the 
wavevectors and group- velocities at the points circled in 
red in the stadium systems presented in Fig. [6j using 
Husimi projections with a coherent spread of Ak/k = 
10%. The states used for the square-lattice system are 
near energies of E = 3.5t, and at E = O.St for the hon- 
eycomb lattice. More details can be found in Appendix 
[A] Fig. [7] shows that the distribution among wavevectors 
emphasizes directions away from the preferred directions 
in group velocity. For certain energy regimes, neither 
wavevectors nor group velocities are evenly distributed 
across all eigenstates of a lattice system. 



D. Internal Bragg Diffraction 

The high-energy eigenstates from Fig. [6] exhibit an un- 
usual behavior: the self-looping classical trajectories that 
are strongly emphasized in the wavefunctions do not ex- 
hibit specular reflection at the boundary. We clarify 
these reflections in the schematics in Fig. [8] Even though 
the absolute angles at each reflection point fall along the 




Figure 8: High-energy states in the square (a) and honeycomb 
(b) lattices can exhibit unusual behaviors, such as group- 
velocity warping and non-specular boundary reflections. The 
former can be seen in the wavefunction (left) by the restric- 
tion of trajectories to 45° diagonals for the square lattice (a) 
and the 60° diagonals for the honeycomb lattice (b). Non- 
specular reflections are magnified in the schematic (right). 
Even though the absolute incoming and outgoing angles for 
each point are the same angle, their angles of incidence (single 
and double arcs) are strongly divergent. 



same diagonal, the angles of incidence vary substantially 
between the incoming and outgoing rays. In the honey- 
comb eigenstate (Fig.^), the reflection consists of scat- 
tering into the other valley and propagating in the exact 
opposite direction. 

While the reflections of many trajectories in high- 
energy states violate specularity as a result of group- 
velocity warping, we have chosen the states in Figs. [6] 
and [8] specifically because these reflections behave in un- 
expected ways. Moreover, these surprising reflections oc- 
cur only at certain points along the boundary where the 
lattice cut deviates from an axis of symmetry; specifically, 
they occur slightly off of clean cuts where jaggedness is 
most prominent. 

To understand these unusual reflections, we use a tech- 
nique called the Gaussian beam[14], which shows the en- 
tire set of wavefunctions available to the system which 
intersect at a particular point in both spatial and mo- 
mentum coordinates. This is accomplished by weighting 
the eigenstates {^e} for a closed system by a coherent 
state |ro,ko,cr) which satisfies the dispersion relation at 
energy Eq. To examine reflections at jagged boundaries, 
we place tq along one of these boundaries and ko point- 




Figure 9: Two Gaussian Beams, constructed by summing the 
set of closed-system eigenstates in the energy range 2.48 < 
E < 2.52 weighted by Eq.[9] using a coherent state with mo- 
mentum uncertainty Ak/k = 5% that sits on the right-hand 
boundary (black circles) with specified momentum (small 
black arrows). The system is a square lattice cut at an 18° 
angle (inset). The incoming group- velocity angle is set to 0° 
at top and —40° at bottom. 



ing away from the bulk. Each eigenstate is associated 
with an eigenenergy E so that the Gaussian beam \I/ is 
defined as 



E 



(^|r ,k ,cr) ij) E . 



(9) 



Because of the finite uncertainty of the coherent state, 
only wavefunctions at energies close to Eo contribute to 
the final result. Thus, only a finite range centered around 
Eq must be considered. 

It is important to choose the spread of the coherent 
state wisely. Too large a coherent state restricts the set of 
eigenstates that contribute to the sum, giving unclear re- 
sults. Too small a coherent state does not provide enough 
information to resolve features of the beam. In Fig. [9] a 
compromise is chosen at Ak/k = 5%, which provides a 
sufficient range of eigenstates to construct a clear beam. 

The classical paths suggested by the Gaussian beams 
in Fig. [9] must all travel through the position ro with 
momentum ftko (Eq. [9| defined by the coherent state 
for each beam. In both top and bottom figures, the 
coherent state lies along the right-hand boundary, al- 
though the wavevectors for each coherent state differs. 
Because the breadth of a coherent state grows in time 
when it propagates, each beam focuses at the coherent 
state, and spread from its center. In both Figs. [9^i and 
|9]3, a specular reflected beam is present, but an addi- 
tional reflected beam emerges as a result internal Bragg 
diffraction [15J. These additional reflections also appear 
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Figure 10: Schematic of internal Bragg diffraction. In light 
gray are the lattice sites for the system, and in dashed outlines 
are shown phantom lattice points outside the system where 
the wavef unction must go to zero. A ray coming in at an 
angle 0i n and reflects at O ut interferes with rays from adjacent 
equivalent points on the boundary unit cell. The angle <j) and 
horizontal distance between adjacent unit cells d are identical 
to the boundary in [9] 



at scattering points for the square lattice in Fig. [8^. 

It is possible to quantify internal Bragg diffraction by 
considering an edge identical to the system in Fig. [9] and 
depicted schematically in Fig. [lOj Here the boundary is 
cut at an angle (j) « 18°, where (j) = is a vertical edge, 



and an incoming plane wave strikes the surface at angle 
#i n , where 0- in = points to the right, and positive angles 
point upward. This plane wave reflects to an outgoing 
angle out where out = points to the left and positive 
angles point upward. 

If there is a repeating unit cell in the edge, two rays 
which hit equivalent points of adjacent boundary unit 
cells gain or lose relative phase based on their wavevec- 
tors and the different distances they travel. For instance, 
a ray incurs an additional phase of S = kd ^^fi when 

#in > —<\> and S = kd S1 " . ° n ^ when in > (j). Here d is the 
horizontal distance between identical points in adjacent 
unit cells and k is the wavevector magnitude of the incom- 
ing wave. When the plane wave is reflected, its neighbor 
gains phase according to the above formulas, but with 
k indicating the outgoing wavevector magnitude. When 
these two phases cancel or add to a multiple of 27r, the 
two rays constructively interfere. Since this is repeated 
over many unit cells, the interference can be quite strong. 

Because the wavelength shrinks with increasing energy, 
more Bragg branches appear as energy goes up. And be- 
cause the distance between adjacent unit cells increases 
for slighter angles against an axis of symmetry, more 
Bragg branches appear for shallower cuts. We find that 
both of these criteria are satisfied for the reflection points 
in Figs. [8] and [5] 

Using the above formulas, we can predict internal 




-3ti/4 -n/2 



n/2 3ti/4 



Si 



in 



Figure 11: The internal Bragg relationship for a square lattice 
with an 18° cut as depicted in Figs. |9| and [To| computed using 
a scattering matrix on a square-lattice with 50 vertical unit 
cells at energy E — 2.5£. The identity line is shown in grey. 
The two incoming group- velocity angles from Fig.[9]of 0° and 
—40° are shown in vertical black dashed lines. The specular 
line is shown in blue, and the upper and lower branches are 
shown in green and red respectively. 



Bragg diffraction for arbitrary cuts and energies. In 
Fig. [TT] we present these results for the system in Fig. [9] 
The boundary unit cell consists of three vertical units 
and one horizontal unit, so that (j) ~ 18°. The two in- 
coming beams from Fig. [9] are represented by vertical 
Each intersects the graph at the 



dashed lines in Fig. 11 



three locations: along the identity line for the incoming 
beam, along the blue specular line for an outgoing beam, 
and along one of the Bragg branches for the other out- 
going beam. Our predictions are strongly validated by 
Fig.0 

The g spread of the test wavepacket used to create the 
Gaussian beam only covers 4 steps along the cut, meaning 
that only a few surface defects can produce substantial 
Bragg scattering. The ubiquity of this effect has impli- 
cations for ray-tracing methods, which bridge classical 
and quantum explanations for phenomena such as fractal 
conductance fluctuations [16| fT7] and caustics[18, 19J and 
encourages a re-examination ray-splitting [20J and other 
hypothetical edge effect s[2T| [22], 

Combining group-velocity restriction and internal 
Bragg diffraction, we argue that the dense linear paths 
in the wavefunctions in Figs. [6] and [8] are indeed linked to 
classical rays which bounce back and forth approximately 
linearly; at one boundary the bounce is non-specular due 
to the cut of the edge and internal Bragg diffraction. 
For the honeycomb lattice, each bounce can be addition- 
ally associated with scattering into another valley. For 
both systems, these wavefunction enhancements are not 
strictly scars[23j, which are generated by unstable classi- 
cal periodic orbits in the analogous classical limit (group 
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velocity) system. Instead, the wavefunction structures 
are likely normal quantum confinement to stable zones 
in classical phase space. 
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IV. CONCLUSIONS 



We have expanded the vector Husimi projection tech- 
nique, introduced in Mason et aZ. [3J , from the continu- 
ous system to lattices. We have demonstrated and ex- 
plained two unusual properties of lattice systems using 
the Husimi projection: group-velocity warping and in- 
ternal Bragg diffraction, both of which can strongly af- 
fect the properties of classical dynamics of these systems, 
in particular producing unexpected self-looping trajecto- 
ries most visible in extreme-energy states. We have also 
shown that Husimi projections can isolate multiple bands 
which are simultaneously represented in the wavefunc- 
tion, using the two valleys of the honeycomb as an ex- 
ample. For the honeycomb lattice, we have shown that 
one can identify locations of scattering between valleys 
by measuring the divergence of the Husimi map for each 
valley separately. 
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Appendix A: Wavevector and Group Velocity 
Distributions 



The plots in Fig. [7] are produced by summing the 
Husimi vectors over many eigenstates of each system in 
Fig. |6j For the square lattice, 600 states with energies 
3A6t < E < 3.54t and for the honeycomb lattice, 300 
states with energies 0.76t < E < 0.84t. This is done for 
256 wavevectors equally separated by angle, and then a 
small Gaussian kernel is applied with angle width 7r/32. 
Each Husimi vector is multiplied by the infinitesimal dk 
determined by the average distance to neighboring vec- 
tors in the sample, and each calculation takes place at 
the points circled in red in Fig. [6] with coherent spread 
of Ak/k = 10%. The contour line in the dispersion re- 
lation is re-computed for each eigenstate to generate the 
coherent states for the Husimi projection. This is done 
to ensure that the steeper gradient of the dispersion re- 
lation near the preferred group velocities does not affect 
our results. 
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